#include "signalCodes.hpp"
#include "topology.hpp"
#include "util.hpp"
#include "meshTypeFunctions.hpp"
#include "mesh_testFunctions.hpp"
#include <catch2/catch_test_macros.hpp>
#include <vector>
#include <utility>

using namespace hms;


TEST_CASE( "Determining non-connecting vertices" )
{
	using hms::otherVertex;
	
	auto mesh { testMeshRectangles() };
	
	SECTION("using direct access via mesh")
	{
		REQUIRE( otherVertex(mesh, 0, 0) == 1 );
		REQUIRE( otherVertex(mesh, 0, 1) == 0 );
		REQUIRE( otherVertex(mesh, 1, 1) == 2 );
		REQUIRE( otherVertex(mesh, 1, 2) == 1 );
		REQUIRE( otherVertex(mesh, 2, 3) == 4 );
		REQUIRE( otherVertex(mesh, 2, 4) == 3 );
		REQUIRE( otherVertex(mesh, 3, 4) == 5 );
		REQUIRE( otherVertex(mesh, 3, 5) == 4 );
		REQUIRE( otherVertex(mesh, 4, 6) == 7 );
		REQUIRE( otherVertex(mesh, 4, 7) == 6 );
		REQUIRE( otherVertex(mesh, 5, 7) == 8 );
		REQUIRE( otherVertex(mesh, 5, 8) == 7 );
		REQUIRE( otherVertex(mesh, 6, 0) == 3 );
		REQUIRE( otherVertex(mesh, 6, 3) == 0 );
		REQUIRE( otherVertex(mesh, 7, 3) == 6 );
		REQUIRE( otherVertex(mesh, 7, 6) == 3 );
		REQUIRE( otherVertex(mesh, 8, 1) == 4 );
		REQUIRE( otherVertex(mesh, 8, 4) == 1 );
		REQUIRE( otherVertex(mesh, 9, 4) == 7 );
		REQUIRE( otherVertex(mesh, 9, 7) == 4 );
		REQUIRE( otherVertex(mesh,10, 2) == 5 );
		REQUIRE( otherVertex(mesh,10, 5) == 2 );
	}
	SECTION("using temporary object hms::Edge")
	{
		REQUIRE( otherVertex(mesh.edge( 0),0) == 1 );
		REQUIRE( otherVertex(mesh.edge( 0),1) == 0 );
		REQUIRE( otherVertex(mesh.edge( 1),1) == 2 );
		REQUIRE( otherVertex(mesh.edge( 1),2) == 1 );
		REQUIRE( otherVertex(mesh.edge( 2),3) == 4 );
		REQUIRE( otherVertex(mesh.edge( 2),4) == 3 );
		REQUIRE( otherVertex(mesh.edge( 3),4) == 5 );
		REQUIRE( otherVertex(mesh.edge( 3),5) == 4 );
		REQUIRE( otherVertex(mesh.edge( 4),6) == 7 );
		REQUIRE( otherVertex(mesh.edge( 4),7) == 6 );
		REQUIRE( otherVertex(mesh.edge( 5),7) == 8 );
		REQUIRE( otherVertex(mesh.edge( 5),8) == 7 );
		REQUIRE( otherVertex(mesh.edge( 6),0) == 3 );
		REQUIRE( otherVertex(mesh.edge( 6),3) == 0 );
		REQUIRE( otherVertex(mesh.edge( 7),3) == 6 );
		REQUIRE( otherVertex(mesh.edge( 7),6) == 3 );
		REQUIRE( otherVertex(mesh.edge( 8),1) == 4 );
		REQUIRE( otherVertex(mesh.edge( 8),4) == 1 );
		REQUIRE( otherVertex(mesh.edge( 9),4) == 7 );
		REQUIRE( otherVertex(mesh.edge( 9),7) == 4 );
		REQUIRE( otherVertex(mesh.edge(10),2) == 5 );
		REQUIRE( otherVertex(mesh.edge(10),5) == 2 );
	}
}


TEST_CASE( "Determining shared vertices" )
{
	using hms::sharedVertex;
	
	auto mesh { testMeshRectangles() };
	
	SECTION( "direct access via mesh and indices" )
	{
		REQUIRE( sharedVertex(mesh, 0,  1) == 1 );
		REQUIRE( sharedVertex(mesh, 0,  6) == 0 );
		REQUIRE( sharedVertex(mesh, 0,  8) == 1 );
		REQUIRE( sharedVertex(mesh, 0,  2) ==
			hms::signalCode::not_found );
		REQUIRE( sharedVertex(mesh, 1,  0) == 1 );
		REQUIRE( sharedVertex(mesh, 1,  8) == 1 );
		REQUIRE( sharedVertex(mesh, 1, 10) == 2 );
		REQUIRE( sharedVertex(mesh, 2,  3) == 4 );
		REQUIRE( sharedVertex(mesh, 2,  6) == 3 );
		REQUIRE( sharedVertex(mesh, 2,  7) == 3 );
		REQUIRE( sharedVertex(mesh, 2,  8) == 4 );
		REQUIRE( sharedVertex(mesh, 2,  9) == 4 );
		REQUIRE( sharedVertex(mesh, 3,  2) == 4 );
		REQUIRE( sharedVertex(mesh, 3,  8) == 4 );
		REQUIRE( sharedVertex(mesh, 3,  9) == 4 );
		REQUIRE( sharedVertex(mesh, 3, 10) == 5 );
		REQUIRE( sharedVertex(mesh, 3, 11) == 5 );
		REQUIRE( sharedVertex(mesh, 4,  5) == 7 );
		REQUIRE( sharedVertex(mesh, 4,  7) == 6 );
		REQUIRE( sharedVertex(mesh, 4,  9) == 7 );
		REQUIRE( sharedVertex(mesh, 5,  4) == 7 );
		REQUIRE( sharedVertex(mesh, 5,  9) == 7 );
		REQUIRE( sharedVertex(mesh, 5, 11) == 8 );
		REQUIRE( sharedVertex(mesh, 6,  0) == 0 );
		REQUIRE( sharedVertex(mesh, 6,  2) == 3 );
		REQUIRE( sharedVertex(mesh, 7,  2) == 3 );
		REQUIRE( sharedVertex(mesh, 7,  4) == 6 );
		REQUIRE( sharedVertex(mesh, 8,  0) == 1 );
		REQUIRE( sharedVertex(mesh, 8,  1) == 1 );
		REQUIRE( sharedVertex(mesh, 8,  2) == 4 );
		REQUIRE( sharedVertex(mesh, 8,  3) == 4 );
		REQUIRE( sharedVertex(mesh, 8,  9) == 4 );
		REQUIRE( sharedVertex(mesh, 9,  2) == 4 );
		REQUIRE( sharedVertex(mesh, 9,  3) == 4 );
		REQUIRE( sharedVertex(mesh, 9,  4) == 7 );
		REQUIRE( sharedVertex(mesh, 9,  5) == 7 );
		REQUIRE( sharedVertex(mesh, 9,  8) == 4 );
		REQUIRE( sharedVertex(mesh,10,  1) == 2 );
		REQUIRE( sharedVertex(mesh,10,  3) == 5 );
		REQUIRE( sharedVertex(mesh,10, 11) == 5 );
		REQUIRE( sharedVertex(mesh,11,  3) == 5 );
		REQUIRE( sharedVertex(mesh,11,  5) == 8 );
		REQUIRE( sharedVertex(mesh,11, 10) == 5 );
	}
	SECTION( "using temporary hms::Edge objects" )
	{
		REQUIRE( sharedVertex(mesh.edge( 0), mesh.edge( 1)) == 1 );
		REQUIRE( sharedVertex(mesh.edge( 0), mesh.edge( 6)) == 0 );
		REQUIRE( sharedVertex(mesh.edge( 0), mesh.edge( 8)) == 1 );
		REQUIRE( sharedVertex(mesh.edge( 0), mesh.edge( 2)) ==
			hms::signalCode::not_found );
		REQUIRE( sharedVertex(mesh.edge( 1), mesh.edge( 0)) == 1 );
		REQUIRE( sharedVertex(mesh.edge( 1), mesh.edge( 8)) == 1 );
		REQUIRE( sharedVertex(mesh.edge( 1), mesh.edge(10)) == 2 );
		REQUIRE( sharedVertex(mesh.edge( 2), mesh.edge( 3)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 2), mesh.edge( 6)) == 3 );
		REQUIRE( sharedVertex(mesh.edge( 2), mesh.edge( 7)) == 3 );
		REQUIRE( sharedVertex(mesh.edge( 2), mesh.edge( 8)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 2), mesh.edge( 9)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 3), mesh.edge( 2)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 3), mesh.edge( 8)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 3), mesh.edge( 9)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 3), mesh.edge(10)) == 5 );
		REQUIRE( sharedVertex(mesh.edge( 3), mesh.edge(11)) == 5 );
		REQUIRE( sharedVertex(mesh.edge( 4), mesh.edge( 5)) == 7 );
		REQUIRE( sharedVertex(mesh.edge( 4), mesh.edge( 7)) == 6 );
		REQUIRE( sharedVertex(mesh.edge( 4), mesh.edge( 9)) == 7 );
		REQUIRE( sharedVertex(mesh.edge( 5), mesh.edge( 4)) == 7 );
		REQUIRE( sharedVertex(mesh.edge( 5), mesh.edge( 9)) == 7 );
		REQUIRE( sharedVertex(mesh.edge( 5), mesh.edge(11)) == 8 );
		REQUIRE( sharedVertex(mesh.edge( 6), mesh.edge( 0)) == 0 );
		REQUIRE( sharedVertex(mesh.edge( 6), mesh.edge( 2)) == 3 );
		REQUIRE( sharedVertex(mesh.edge( 7), mesh.edge( 2)) == 3 );
		REQUIRE( sharedVertex(mesh.edge( 7), mesh.edge( 4)) == 6 );
		REQUIRE( sharedVertex(mesh.edge( 8), mesh.edge( 0)) == 1 );
		REQUIRE( sharedVertex(mesh.edge( 8), mesh.edge( 1)) == 1 );
		REQUIRE( sharedVertex(mesh.edge( 8), mesh.edge( 2)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 8), mesh.edge( 3)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 8), mesh.edge( 9)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 9), mesh.edge( 2)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 9), mesh.edge( 3)) == 4 );
		REQUIRE( sharedVertex(mesh.edge( 9), mesh.edge( 4)) == 7 );
		REQUIRE( sharedVertex(mesh.edge( 9), mesh.edge( 5)) == 7 );
		REQUIRE( sharedVertex(mesh.edge( 9), mesh.edge( 8)) == 4 );
		REQUIRE( sharedVertex(mesh.edge(10), mesh.edge( 1)) == 2 );
		REQUIRE( sharedVertex(mesh.edge(10), mesh.edge( 3)) == 5 );
		REQUIRE( sharedVertex(mesh.edge(10), mesh.edge(11)) == 5 );
		REQUIRE( sharedVertex(mesh.edge(11), mesh.edge( 3)) == 5 );
		REQUIRE( sharedVertex(mesh.edge(11), mesh.edge( 5)) == 8 );
		REQUIRE( sharedVertex(mesh.edge(11), mesh.edge(10)) == 5 );
	}
}


TEST_CASE( "Finding edges attached to each vertex" )
{
	SECTION("rectangular mesh")
	{
		auto mesh { testMeshRectangles() };
		
		auto edgesAtVertex { hms::edgesAtVertices( mesh ) };
		
		REQUIRE( edgesAtVertex[0] == IdxVector { 0,  6        } );
		REQUIRE( edgesAtVertex[1] == IdxVector { 0,  1,  8    } );
		REQUIRE( edgesAtVertex[2] == IdxVector { 1, 10        } );
		REQUIRE( edgesAtVertex[3] == IdxVector { 2,  6,  7    } );
		REQUIRE( edgesAtVertex[4] == IdxVector { 2,  3,  8, 9 } );
		REQUIRE( edgesAtVertex[5] == IdxVector { 3, 10, 11    } );
		REQUIRE( edgesAtVertex[6] == IdxVector { 4,  7        } );
		REQUIRE( edgesAtVertex[7] == IdxVector { 4,  5,  9    } );
		REQUIRE( edgesAtVertex[8] == IdxVector { 5, 11        } );
	}
	SECTION("triangular mesh")
	{
		auto mesh { testMeshTriangles() };
		
		auto edgesAtVertex { hms::edgesAtVertices( mesh ) };
		
		REQUIRE( edgesAtVertex[0] == IdxVector { 0, 1, 2, 3, 4 } );
		REQUIRE( edgesAtVertex[1] == IdxVector { 4, 5          } );
		REQUIRE( edgesAtVertex[2] == IdxVector { 3, 5, 6       } );
		REQUIRE( edgesAtVertex[3] == IdxVector { 2, 6, 7       } );
		REQUIRE( edgesAtVertex[4] == IdxVector { 1, 7, 8       } );
		REQUIRE( edgesAtVertex[5] == IdxVector { 0, 8          } );
	}
	SECTION("mixed mesh")
	{
		auto mesh { testMeshMixed() };
		
		auto edgesAtVertex { hms::edgesAtVertices( mesh ) };
		
		REQUIRE( edgesAtVertex[0] == IdxVector { 0, 1    } );
		REQUIRE( edgesAtVertex[1] == IdxVector { 0, 2, 3 } );
		REQUIRE( edgesAtVertex[2] == IdxVector { 1, 2, 4 } );
		REQUIRE( edgesAtVertex[3] == IdxVector { 3, 5    } );
		REQUIRE( edgesAtVertex[4] == IdxVector { 4, 5    } );
	}
}


TEST_CASE( "Finding closed loops formed by edges, rectangular mesh" )
{
	using hms::findLoops;
	using hms::sameEntries;
	using hms::edgesAtVertices;
	
	SECTION("rectangular mesh, access via indices")
	{
		auto mesh { testMeshRectangles() };
		
		auto edgesAtVertex { edgesAtVertices( mesh ) };
		
		REQUIRE( sameEntries(
			findLoops(mesh, 0, edgesAtVertex)[0],
			IdxVector{ 0, 6, 2, 8} ) );
		REQUIRE( sameEntries(
			findLoops(mesh, 1,edgesAtVertex)[0],
			IdxVector{ 1, 8, 3,10} ) );
		REQUIRE( sameEntries(
			findLoops(mesh, 2,edgesAtVertex)[0],
			IdxVector{ 2, 7, 4, 9} ) );
		REQUIRE( sameEntries(
			findLoops(mesh, 3,edgesAtVertex)[0],
			IdxVector{ 3, 9, 5,11} ) );
		REQUIRE( findLoops(mesh,  4, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,  5, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,  6, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,  7, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,  8, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,  9, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 10, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 11, edgesAtVertex).empty() );
	}
	SECTION("rectangular mesh, using temporary hms::Edge object")
	{
		auto mesh { testMeshRectangles() };
		
		auto edgesAtVertex { edgesAtVertices( mesh ) };
		
		REQUIRE( sameEntries(
			findLoops(mesh, mesh.edge(0),edgesAtVertex)[0],
			IdxVector{ 0, 6, 2, 8} ) );
		REQUIRE( sameEntries(
			findLoops(mesh, mesh.edge(1),edgesAtVertex)[0],
			IdxVector{ 1, 8, 3,10} ) );
		REQUIRE( sameEntries(
			findLoops(mesh, mesh.edge(2),edgesAtVertex)[0],
			IdxVector{ 2, 7, 4, 9} ) );
		REQUIRE( sameEntries(
			findLoops(mesh, mesh.edge(3),edgesAtVertex)[0],
			IdxVector{ 3, 9, 5,11} ) );
		REQUIRE( findLoops(mesh,mesh.edge( 4),edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,mesh.edge( 5),edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,mesh.edge( 6),edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,mesh.edge( 7),edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,mesh.edge( 8),edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,mesh.edge( 9),edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,mesh.edge(10),edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh,mesh.edge(11),edgesAtVertex).empty() );
	}
	SECTION("triangular mesh, direct access via indices")
	{
		auto mesh { testMeshTriangles() };
		
		auto edgesAtVertex { edgesAtVertices( mesh ) };
		
		REQUIRE( sameEntries(
			findLoops(mesh, 0, edgesAtVertex)[0],
			IdxVector{ 0, 1, 8 } ) );
		REQUIRE( sameEntries(
			findLoops(mesh, 1, edgesAtVertex)[0],
			IdxVector{ 1, 2, 7 } ) );
		REQUIRE( sameEntries(
			findLoops(mesh, 2, edgesAtVertex)[0],
			IdxVector{ 2, 3, 6 } ) );
		REQUIRE( sameEntries(
			findLoops(mesh, 3, edgesAtVertex)[0],
			IdxVector{ 3, 4, 5 } ) );
		REQUIRE( findLoops(mesh, 4, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 5, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 6, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 7, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 8, edgesAtVertex).empty() );
	}
	SECTION( "mixed mesh, direct access" )
	{
		auto mesh { testMeshMixed() };
		
		auto edgesAtVertex { edgesAtVertices( mesh ) };
		
		REQUIRE( sameEntries(
			findLoops(mesh, 0, edgesAtVertex)[0],
			IdxVector{ 0, 1, 2 } ));
		REQUIRE( findLoops(mesh, 1, edgesAtVertex).empty() );
		REQUIRE( sameEntries(
			findLoops(mesh, 2, edgesAtVertex)[0],
			IdxVector{ 2, 3, 4, 5 } ) );
		REQUIRE( findLoops(mesh, 3, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 4, edgesAtVertex).empty() );
		REQUIRE( findLoops(mesh, 5, edgesAtVertex).empty() );
	}
}


TEST_CASE( "edgeLoopToVertices: Retrieving unique vertices from edge loop")
{
	using hms::sameEntries;
	using hms::edgeLoopToVertices;
	
	
	SECTION( "rectangular mesh" )
	{
		auto mesh { testMeshRectangles(true) };
		
		Array4i edgeLoop;
		edgeLoop = { 6, 2, 8, 0 };
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, edgeLoop ),
			Array4i { 0, 3, 4, 1 } ));
		
		edgeLoop << 0, 8, 2, 6;
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, edgeLoop ),
			IdxVector { 0, 1, 4, 3 } ));
		
		edgeLoop << 0, 2, 6, 8;
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, edgeLoop ),
			IdxVector { 0, 1, 4, 3 } ));
	}
	SECTION( "rectangular mesh, passing cell index" )
	{
		auto mesh { testMeshRectangles(true) };
		
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, 0 ),
			Array4i { 0, 3, 4, 1 } ));
	}
	
	SECTION( "triangular mesh" )
	{
		auto mesh { testMeshTriangles(true) };
		
		Array3i edgeLoop;
		edgeLoop = { 3, 8, 1 };
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, edgeLoop ),
			Array3i { 0, 5, 4 } ));
		
		edgeLoop = { 3, 1, 8 };
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, edgeLoop ),
			Array3i { 5, 0, 4 } ));
	}
	SECTION( "triangular mesh, passing cell index" )
	{
		auto mesh { testMeshTriangles(true) };

		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, 0 ),
			Array3i { 0, 5, 4 } ));
	}
	SECTION( "mixed mesh" )
	{
		auto mesh { testMeshMixed(true) };
		
		ArrayXi edgeLoop (3);
		edgeLoop << 2, 0, 1;
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, edgeLoop ),
			Array3i { 0, 1, 2 } ));
		
		edgeLoop.resize(4);
		edgeLoop << 0, 3, 4, 5;
		REQUIRE( sameEntries(
			edgeLoopToVertices( mesh, edgeLoop ),
			Array4i { 1, 2, 3, 4 } ));
	}
		
}

TEST_CASE( "Determining on which Side of a line a point lies"
){
	using hms::whichSide;
	using hms::Side;
	
	SECTION( "three vectors" )
	{
		std::vector<Vector2s> line { {1,1}, {2,8} };
		
		REQUIRE( whichSide( line[0], line[1], Vector2s{1,2} ) == Side::left );
		REQUIRE( whichSide( line[0], line[1], Vector2s{1,0} ) == Side::right );
		REQUIRE( whichSide( line[0], line[1], Vector2s{0,-6} ) == Side::on );
		
		REQUIRE( whichSide( line[1], line[0], Vector2s{1,2} ) == Side::right );
		REQUIRE( whichSide( line[1], line[0], Vector2s{1,0} ) == Side::left );
		REQUIRE( whichSide( line[1], line[0], Vector2s{0,-6} ) == Side::on );
	}
	SECTION( "Mesh, edge index and vector" )
	{
		hms::Mesh mesh { testMeshMixed() };
		
		REQUIRE( whichSide( mesh, 0, Vector2s{0,0} ) == Side::left  );
		REQUIRE( whichSide( mesh, 1, Vector2s{0,0} ) == Side::right );
		REQUIRE( whichSide( mesh, 2, Vector2s{0,0} ) == Side::right );
		REQUIRE( whichSide( mesh, 3, Vector2s{0,0} ) == Side::on    );
		REQUIRE( whichSide( mesh, 4, Vector2s{0,0} ) == Side::right );
		REQUIRE( whichSide( mesh, 5, Vector2s{0,0} ) == Side::on    );
		
		REQUIRE( whichSide( mesh, 0, Vector2s{0.5,0.5} ) == Side::left  );
		REQUIRE( whichSide( mesh, 1, Vector2s{0.5,0.5} ) == Side::right );
		REQUIRE( whichSide( mesh, 2, Vector2s{0.5,0.5} ) == Side::right );
		REQUIRE( whichSide( mesh, 3, Vector2s{0.5,0.5} ) == Side::left  );
		REQUIRE( whichSide( mesh, 4, Vector2s{0.5,0.5} ) == Side::right );
		REQUIRE( whichSide( mesh, 5, Vector2s{0.5,0.5} ) == Side::left  );
	}
	SECTION( "Edge and vector" )
	{
		hms::Mesh mesh { testMeshMixed() };
		
		REQUIRE( whichSide( mesh.edge(0), Vector2s{0,0} ) == Side::left  );
		REQUIRE( whichSide( mesh.edge(1), Vector2s{0,0} ) == Side::right );
		REQUIRE( whichSide( mesh.edge(2), Vector2s{0,0} ) == Side::right );
		REQUIRE( whichSide( mesh.edge(3), Vector2s{0,0} ) == Side::on    );
		REQUIRE( whichSide( mesh.edge(4), Vector2s{0,0} ) == Side::right );
		REQUIRE( whichSide( mesh.edge(5), Vector2s{0,0} ) == Side::on    );
		
		REQUIRE( whichSide( mesh.edge(0), Vector2s{0.5,0.5} ) == Side::left  );
		REQUIRE( whichSide( mesh.edge(1), Vector2s{0.5,0.5} ) == Side::right );
		REQUIRE( whichSide( mesh.edge(2), Vector2s{0.5,0.5} ) == Side::right );
		REQUIRE( whichSide( mesh.edge(3), Vector2s{0.5,0.5} ) == Side::left  );
		REQUIRE( whichSide( mesh.edge(4), Vector2s{0.5,0.5} ) == Side::right );
		REQUIRE( whichSide( mesh.edge(5), Vector2s{0.5,0.5} ) == Side::left  );
	}
}

TEST_CASE( "Bottom left index is searched" )
{
	using hms::findBottomLeft;
	
	Matrix2Xs verts1 (2,4);
	verts1 <<
		1, 0, 1, 0,
		0, 0, 1, 1;
	
	Matrix2Xs verts2 (2,4);
	verts2 <<
		-20.001, -20    ,  20    , 20,
		 20    , -20.001, -20.002, 20;
	
	Matrix2Xs verts3 (2,3);
	verts3 <<
		-1, 1, 0,
		-1,-1, 1;
	
	hms::Mesh mesh1 { std::move(verts1) };
	hms::Mesh mesh2 { std::move(verts2) };
	hms::Mesh mesh3 { std::move(verts3) };
	
	SECTION( "Quadrilaterals" )
	{
		REQUIRE( findBottomLeft( mesh1, Array4i{0,1,2,3} ) == 1 );
		REQUIRE( findBottomLeft( mesh2, Array4i{0,1,2,3} ) == 2 );
	}
	SECTION( "Triangles" )
	{	
		REQUIRE( findBottomLeft( mesh3, Array3i{0,1,2} ) == 0 );
	}
}


TEST_CASE( "Ordering vertices" )
{
	using hms::orderVertices;
	
	{
	Matrix2Xs verts (2, 4);
	verts <<
		-3,  1, -2, 3,
		-1, -4,  3, 1;
	
	Matrix2Xs vertsCopy { verts.replicate(1,1) };
	hms::Mesh mesh { std::move( vertsCopy ) };
	SECTION( "single quadrilateral" )
	{
		Array4i vertInds {0,1,2,3};
		
		/* order vertices clockwise */
		orderVertices( mesh, vertInds, hms::VertexOrdering::clockwise );
		
		/* check indices */
		REQUIRE( (vertInds == Array4i{1,0,2,3}).all() );
		
		/* check vertices behind indices */
		for (int i {0}; i<4; ++i)
			REQUIRE( mesh.vertex(vertInds[i]) == verts.col(vertInds[i]) );
		
		/* order vertices counterclockwise */
		orderVertices(mesh, vertInds, hms::VertexOrdering::counterclockwise);
		REQUIRE( (vertInds == Array4i{1,3,2,0}).all() );
	}
	}
	SECTION( "standard test mesh (rect)" )
	{
		hms::Mesh mesh { testMeshRectangles() };
		Array4i verts;
		Array4i cwise;
		Array4i ccwise;
		
		/* bottom left cell */
		verts  << 4, 3, 0, 1;
		cwise  << 0, 1, 4, 3;
		ccwise << 0, 3, 4, 1;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* bottom right */
		verts  << 7, 6, 3, 4;
		cwise  << 3, 4, 7, 6;
		ccwise << 3, 6, 7, 4;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* top left */
		verts  << 5, 1, 4, 2;
		cwise  << 1, 2, 5, 4;
		ccwise << 1, 4, 5, 2;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* top right cell */
		verts  << 5, 8, 7, 4;
		cwise  << 4, 5, 8, 7;
		ccwise << 4, 7, 8, 5;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* top right cell, in order */
		verts  << 4, 5, 8, 7;
		cwise  << 4, 5, 8, 7;
		ccwise << 4, 7, 8, 5;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* outer corners */
		verts  << 6, 8, 2, 0;
		cwise  << 0, 2, 8, 6;
		ccwise << 0, 6, 8, 2;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
	}
	SECTION( "triangular cells" )
	{
		auto mesh { testMeshTriangles() };
		Array3i verts;
		Array3i cwise;
		Array3i ccwise;
		
		/* cell 0 */
		verts  << 5, 4, 0;
		cwise  << 0, 4, 5;
		ccwise << 0, 5, 4;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* cell 1 */
		verts  << 3, 4, 0;
		cwise  << 0, 3, 4;
		ccwise << 0, 4, 3;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* cell 2 */
		verts  << 2, 3, 0;
		cwise  << 0, 2, 3;
		ccwise << 0, 3, 2;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
		
		/* cell 3 */
		verts  << 0, 1, 2;
		cwise  << 1, 2, 0;
		ccwise << 1, 0, 2;
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		REQUIRE( (verts == cwise).all() );
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		REQUIRE( (verts == ccwise).all() );
	}
}


TEST_CASE( "Ordering a loop of edges" )
{
	using hms::orderVertices;
	using hms::orderEdgeLoop;
	using hms::edgeLoopToVertices;
	
	ArrayXi edgeLoop;
	ArrayXi orderedLoop;
	ArrayXi verts;
	
	SECTION( "rectangular mesh" )
	{
		auto mesh { testMeshRectangles(true) };
		edgeLoop.resize(4);
		orderedLoop.resize(4);
		
		edgeLoop << 6, 2, 0, 8;
		orderedLoop = edgeLoop;
		verts = edgeLoopToVertices(mesh, edgeLoop);
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		orderEdgeLoop(mesh, verts, edgeLoop);
		REQUIRE( (edgeLoop == orderedLoop).all() );
		
		edgeLoop << 0, 8, 2, 6 ;
		verts = edgeLoopToVertices(mesh, edgeLoop);
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		orderEdgeLoop(mesh, verts, edgeLoop);
		REQUIRE( (edgeLoop == orderedLoop).all() );
		
		edgeLoop << 0, 2, 6, 8 ;
		orderedLoop << 8, 0, 2, 6 ;
		verts = edgeLoopToVertices(mesh, edgeLoop);
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		orderEdgeLoop(mesh, verts, edgeLoop);
		REQUIRE( (edgeLoop == orderedLoop).all() );
		
		edgeLoop << 2, 7, 4, 1 ;
		orderedLoop << 2, 1, 4, 7 ;
		verts = edgeLoopToVertices(mesh, edgeLoop);
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		orderEdgeLoop(mesh, verts, edgeLoop);
		REQUIRE( (edgeLoop == orderedLoop).all() );
	}
	SECTION( "triangular mesh" )
	{
		auto mesh { testMeshTriangles(true) };
		edgeLoop.resize(3);
		orderedLoop.resize(3);
		
		edgeLoop << 4, 0, 5 ;
		orderedLoop = edgeLoop;
		verts = edgeLoopToVertices(mesh, edgeLoop);
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		orderEdgeLoop(mesh, verts, edgeLoop);
		REQUIRE( (edgeLoop == orderedLoop).all() );
		
		edgeLoop << 5, 4, 0 ;
		verts = edgeLoopToVertices(mesh, edgeLoop);
		orderVertices(mesh, verts, hms::VertexOrdering::counterclockwise);
		orderEdgeLoop(mesh, verts, edgeLoop);
		REQUIRE( (edgeLoop == orderedLoop).all() );
		
		edgeLoop << 7, 1, 2 ;
		orderedLoop << 2, 7, 1 ;
		verts = edgeLoopToVertices(mesh, edgeLoop);
		orderVertices(mesh, verts, hms::VertexOrdering::clockwise);
		orderEdgeLoop(mesh, verts, edgeLoop);
		REQUIRE( (edgeLoop == orderedLoop).all() );
	}
}

TEST_CASE( "connectingEdges: Finding edges which connect given vertices" ){
	hms::Mesh mesh { testMeshRectangles(true) };

	std::vector<IdxVector> edgesAtVertex { hms::edgesAtVertices(mesh) };
	
	using hms::connectingEdges;
	REQUIRE( connectingEdges(mesh, {0,3,6}, edgesAtVertex) == IdxVector{6,7} );
	REQUIRE( connectingEdges(mesh, {1,4,5,8}, edgesAtVertex) == IdxVector{0,3,11});
}

TEST_CASE( "Determining bounding box of polygon" )
{
	using hms::boundingBox;
	
	Vector2s a { 0  , 0.5};
	Vector2s b { 0.7, 0  };
	Vector2s c { 1  , 0.2};
	Vector2s d { 0.9, 0.8};
	Vector2s e { 0.2, 1  };
	
	SECTION( "Using polygon corners stored as column vectors in matrix" )
	{
		Matrix2Xs polygon ( hms::dims, 5 );
		polygon << a, b, c, d, e;
		
		AlignedBox box { boundingBox(polygon) };
		REQUIRE( box.corner(AlignedBox::BottomLeft) == Vector2s{0,0} );
		REQUIRE( box.corner(AlignedBox::TopRight  ) == Vector2s{1,1} );
	}
	SECTION( "Using polygon corners stored in std::vector" )
	{
		std::vector<Vector2s> polygon { a, b, c, d, e };
		
		AlignedBox box { boundingBox(polygon) };
		REQUIRE( box.corner(AlignedBox::BottomLeft) == Vector2s{0,0} );
		REQUIRE( box.corner(AlignedBox::TopRight  ) == Vector2s{1,1} );
	}
}

TEST_CASE( "Testing behaviour of AlignedBox::contains" )
{	
	Vector2s a {0,0}, b {1,1};
	AlignedBox bb { a, b };
	
	REQUIRE( bb.contains(a) );
	REQUIRE( bb.contains(b) );
	REQUIRE( bb.contains( Vector2s{ 0.5, 0.5} ) == true );
	REQUIRE( bb.contains( Vector2s{ 1.1, 0  } ) == false);
	REQUIRE( bb.contains( Vector2s{ 0  , 1.1} ) == false);
	REQUIRE( bb.contains( Vector2s{-0.1, 0  } ) == false);
	REQUIRE( bb.contains( Vector2s{ 0  ,-0.1} ) == false);
}

TEST_CASE( "Determining whether point lies in polygon" )
{
	using hms::inPolygon;
	
	/* self-intersecting polygon which can cover most test scenarios */
	std::vector<Vector2s> polygon {
		{ 0, 0},
		{ 4, 4},
		{ 0, 4},
		{ 5,-1},
		{ 8, 5},
		{ 7, 8},
		{ 5, 8},
		{ 3, 5},
		{ 1,10},
		{-3, 5}
	};
	
	std::vector<Vector2s> points {
		{-1, 4},
		{ 1, 8},
		{ 2, 3},
		{ 2, 4},
		{ 2, 5},
		{ 6, 5},
		{ 2, 1},
		{ 3, 6},
		{ 7, 0},
		{-2, 0}
	};
	
	REQUIRE( inPolygon( polygon, points[0] ) == true );
	REQUIRE( inPolygon( polygon, points[1] ) == true );
	REQUIRE( inPolygon( polygon, points[2] ) == true );
	REQUIRE( inPolygon( polygon, points[3] ) == true );
	REQUIRE( inPolygon( polygon, points[4] ) == true );
	REQUIRE( inPolygon( polygon, points[5] ) == true );
	REQUIRE( inPolygon( polygon, points[6] ) == false);
	REQUIRE( inPolygon( polygon, points[7] ) == false);
	REQUIRE( inPolygon( polygon, points[8] ) == false);
	REQUIRE( inPolygon( polygon, points[9] ) == false);
}
		

TEST_CASE( "Finding points which lie inside a polygon" )
{
	using hms::sameEntries;
	using hms::verticesInPolygon;
	
	auto mesh { testMeshRectangles() };
	
	std::vector<Vector2s> polygon { {0,2}, {1,2}, {0,3} };
	REQUIRE( verticesInPolygon( mesh, polygon ).empty() );
	
	/* triangle, check points inside */
	polygon = { {-1.1, -1.1}, {1.1, 0}, {0, 0.1} };
	REQUIRE( verticesInPolygon( mesh, polygon ) == IdxVector{ 0,4,7 } );
			
	/* same triangle, sorted differently */
	polygon = { {-1.1, -1.1}, {0, 0.1}, {1.1, 0} };
	REQUIRE( verticesInPolygon( mesh, polygon ) == IdxVector{ 0,4,7 } );
	
	/* concave polygon around top and right boundary */
	polygon = { {1, -1.1}, {1.1, 1.1}, {-1.1, 1}, {0.9, 0.9} };
	REQUIRE( verticesInPolygon( mesh, polygon ) == IdxVector{ 2,5,6,7,8 } );
}

TEST_CASE( "Finding edges which lie inside a polygon" )
{
	using hms::sameEntries;
	using hms::edgesInPolygon;
	
	auto mesh { testMeshRectangles() };
	
	std::vector<Vector2s> polygon { {0,2}, {1,2}, {0,3} };
	REQUIRE( edgesInPolygon( mesh, polygon ).empty() );
	
	/* concave polygon around top and right boundary */
	polygon = { {1, -1.1}, {1.1, 1.1}, {-1.1, 1}, {0.9, 0.9} };
	REQUIRE( sameEntries(
		edgesInPolygon( mesh, polygon ), IdxVector{ 4,5,10,11 } ) );
}

TEST_CASE( "Finding cells whose centroids lie inside a polygon" ){
	using hms::sameEntries;
	using hms::cellsInPolygon;
	
	auto mesh { testMeshRectangles(true) };
	
	/* polygon looks something like this, with a self-intersection at (0,0): 
	 *  /|
	 * /_|__
	 *   | /
	 *   |/
	 */
	std::vector<Vector2s> polygon{
		{ 0,-2},
		{ 1, 0},
		{-1, 0},
		{ 0, 2} };
	
	REQUIRE( cellsInPolygon( mesh, polygon ) == IdxVector{1,2} );
}


TEST_CASE( "Finding vertices which lie inside a rectangle" )
{
	using hms::sameEntries;
	using hms::verticesInBox;
	
	auto mesh { testMeshRectangles(true) };
	
	REQUIRE( sameEntries(
		verticesInBox(mesh, {-1.1,-1.1}, {-0.9, 1.1} ),
		IdxVector{0,1,2} ) );
	
	REQUIRE( sameEntries(
		verticesInBox(mesh, {-1.1,-1.1}, { 1.1, 1.1} ),
		IdxVector{0,1,2,3,4,5,6,7,8} ) );
	
	REQUIRE( sameEntries(
		verticesInBox(mesh, {-0.1, 0.9}, { 1.1, 1.1} ),
		IdxVector{5,8} ) );
	
	REQUIRE( sameEntries(
		verticesInBox(mesh, {-0.1,-0.1}, { 1.1, 1.1} ),
		IdxVector{4,5,7,8} ) );
	
	REQUIRE( sameEntries(
		verticesInBox(mesh, { 1.1, 1.1}, {-0.1,-0.1} ),
		IdxVector{4,5,7,8} ) );
	
	REQUIRE( sameEntries(
		verticesInBox(mesh, { 0.9, 1.1}, { 1.1, 0.9} ),
		IdxVector{8} ) );
	
	REQUIRE( sameEntries(
		verticesInBox(mesh, { 0,-1}, { 0, 1} ),
		IdxVector{3,4,5} ) );
}

TEST_CASE( "Finding edges which lie inside a rectangle" )
{
	using hms::sameEntries;
	using hms::edgesInBox;
	
	auto mesh { testMeshRectangles(true) };
	
	REQUIRE( sameEntries(
		edgesInBox(mesh, {-1.1,-1.1}, {-0.9, 1.1} ),
		IdxVector{8,9} ) );
	
	REQUIRE( sameEntries(
		edgesInBox(mesh, {-1.1,-1.1}, { 1.1, 1.1} ),
		IdxVector{0,1,2,3,4,5,6,7,8,9,10,11} ) );
	
	REQUIRE( edgesInBox(mesh, {-0.1, 0.9}, { 1.1, 1.1} ) == IdxVector{11} );
	
	REQUIRE( sameEntries(
		edgesInBox(mesh, {-0.1,-0.1}, { 1.1, 1.1} ),
		IdxVector{1,3,5,11} ) );
	
	REQUIRE( sameEntries(
		edgesInBox(mesh, { 1.1, 1.1}, {-0.1,-0.1} ),
		IdxVector{1,3,5,11} ) );
	
	REQUIRE( edgesInBox(mesh, { 0.9, 1.1}, { 1.1, 0.9} ).empty() );
	
	REQUIRE( sameEntries(
		edgesInBox(mesh, { 0,-1}, { 0, 1} ),
		IdxVector{2,3} ) );
}

TEST_CASE( "Finding cells which lie inside a rectangle" )
{
	using hms::cellsInBox;
	
	auto mesh { testMeshTriangles(true) };
	
	REQUIRE( cellsInBox(mesh, {1.1,2.1}, {-1.4,0} ) == IdxVector{0,1,2} );
}



TEST_CASE( "Tracing line of edges through vertex list" )
{
	using hms::sameEntries;
	using hms::verticesToEdgeLine;
	
	auto mesh { testMeshRectangles(true) };
	auto edgesAtVertex { hms::edgesAtVertices(mesh) };
	
	REQUIRE(
		verticesToEdgeLine(mesh, {0,1,2}, edgesAtVertex ) == 
		IdxVector {8,9} );
	REQUIRE(
		verticesToEdgeLine(mesh, {0,1,4}, edgesAtVertex ) ==
		IdxVector {8,0} );
	REQUIRE(
		verticesToEdgeLine(mesh, {7,4,3}, edgesAtVertex ) ==
		IdxVector {1,2} );
	REQUIRE( 
		verticesToEdgeLine(mesh, {6,7,8,5,2}, edgesAtVertex ) == 
		IdxVector {4,5,11,10} );
	REQUIRE( 
		verticesToEdgeLine(mesh, {6,7,4,5,2}, edgesAtVertex ) == 
		IdxVector {4,1,3,10} );
}

TEST_CASE("Finding cell containing point"){

	auto test = [](auto&& dummyMesh){
		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		auto mesh { uniMeshAs<MeshType>(3, 4, 1.5, 0.5, true, -2, 0.5) };
		std::optional<Index> cell {std::nullopt};
		auto reqVal = [&](scalar x, scalar y, Index i, Index j){
			cell = containingCell(mesh, {x, y});
			REQUIRE( cell.has_value() );
			REQUIRE( cell.value() == mesh.cell(i,j).index() );
		};
		auto noVal = [&](scalar x, scalar y){
			cell = containingCell(mesh, {x, y});
			REQUIRE( !cell.has_value() );
		};
		/* inside cell */
		reqVal(-1, 0.75, 0, 0);
		reqVal(1.5, 1.25, 2, 1);
		/* on edge */
		reqVal(-1, 1, 0, 1);
		/* ghost */
		reqVal(0, 2.75, 1, 4);
		/* ghost, on edge */
		reqVal(2.5, 0.75, 3, 0);
		/* on corner */
		reqVal(1, 1, 2, 1);
		/* outside */
		noVal(0, -0.1);
		noVal(-4, 1);
		noVal(0, 3.1);
		noVal(4.1, 1);
	};
	SECTION("UniMesh"){
		test(UniMesh{});
	}
	SECTION("RectMesh"){
		test(RectMesh{});
	}
}

TEST_CASE("Cells along line"){

	auto test = [](auto&& dummyMesh){
		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		auto mesh { uniMeshAs<MeshType>(3, 4, 1.5, 0.5, true, -2, 0.5) };

		REQUIRE( sameEntries(
			cellsAlongLine(mesh, {-0.6, 0.6}, {2, 1.6}), IdxVector{6,7,12,13,18}
		));
		REQUIRE( sameEntries(
			cellsAlongLine(mesh, {2, 1.6}, {-0.6, 0.6}), IdxVector{6,7,12,13,18}
		));
		REQUIRE( sameEntries(
			cellsAlongLine(mesh, {0.5, 2.1}, {0.6, 0.6}), IdxVector{22,17,12,7}
		));
		REQUIRE( sameEntries(
			cellsAlongLine(mesh, {0.6, 0.6}, {0.5, 2.1}), IdxVector{22,17,12,7}
		));
		/* 9 is a ghost cell, which makes the algorithm break off,
		 * so reversing the points does not give the same result */
		REQUIRE( sameEntries(
			cellsAlongLine(mesh, {-1, 0.7}, {3, 0.8}), IdxVector{6,7,8,9}
		));
		REQUIRE( sameEntries(
			cellsAlongLine(mesh, {0.5, 2.1}, {0.6, 0.6}), IdxVector{22,17,12,7}
		));
	};
	SECTION("UniMesh"){
		test(UniMesh{});
	}
	SECTION("RectMesh"){
		test(RectMesh{});
	}
}